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ABSTRACT 

We use a series of high-resolution simulations of a 'Milky- Way' halo coupled to semi- 
analytic methods to study the formation of our own Galaxy and of its stellar halo. 
The physical properties of our model Milky Way, as well as the age and metallicity 
distribution of stars in the different components, are in relatively good agreement with 
observational measurements. Assuming that the stellar halo builds up from the cores 
of the satellite galaxies that merged with the Milky Way over its life-time, we are able 
to study the physical and structural properties of this component. In agreement with 
previous work, we find that the largest contribution to the stellar halo should come 
from a few relatively massive (10* — 10^° Mq) satellites, accreted at early times. Our 
"stellar halo" does not exhibit any clear metallicity gradient, but higher metallicity 
stars are more centrally concentrated than stars of lower abundance. This indicates 
that the probability of observing low-metallicity halo stars increases with distance 
from the Galactic centre. We find that the proposed "dual" nature of the Galactic 
stellar halo can be explained in our model as a result of a mass-metallicity relation 
imprinted in the building blocks of this component. 



Key words: 

Galaxy: halo 



Galaxy: formation - Galaxy: evolution - Galaxy: stellar content 



1 INTRODUCTION 

Our own galaxy - the Milky Way - is a fairly large spiral 
galaxy consisting of four main stellar components. Most of 
the stars are distributed in a thin disk, exhibit a wide range 
of ages, and are on high angular momentum orbits. A much 
smaller mass of stars (about 10-20 per cent of that in the thin 
disk) reside in a distinct component which was established 
abou t 25 years ago though star counts (jGilmore fc ReidI 
119831 ). and which is referred to as the 'thick disk'. The stars 
in the thick disk are old, have on average lower metallicity 
than those of similar age in the thin disk, and are on orbits of 
lower angular momentum. The Galactic bulge is dominated 
by an old and relatively metal-rich stellar popul ation with a 
tail t o low abundances. It has a peanut-shape (jDwek et al.l 
1 19951 ). is kinematically hotter t han the Milky Way dis k but 
colder than the Milky Way halo (JKuiiken fc Richll2002l ). The 
stellar halo represents only a tiny fraction of the total s tellar 
mass (~ 2 X 10^ M© - ICarnev. Latham fc Lairdlll989l ). and 
is dominated by old and metal poor stars which reside on 
low angular momentum orbits. 

While the Milky Way is only one galaxy, it is the one 
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that we can study in unique detail. Accurate measurements 
of ages, metallicities, and kinematics have been collected 
over the years for a large number of individual stars. Over 
the next decade, a number of astrometric and spectroscopic 
surveys will provide accurate spatial, kinematic, and chem- 
ical information for a much larger numb er of stars (e.g. 
the satellite Gaia - IPerrvman et al. 2001; the Radial Ve- 
locity Experiment RAVE - ISteinmetz et al. 2006 : and the 
Sloan Extension for Galactic Und erstanding and Explo- 
ration SEGUE - iBeers et all |2004 ). This vast amount of 
fossil information will provide important advances in our 
understanding of the se quence of events which led to the for- 
mation of our Galaxy (|Freeman fc Bland-Hawthornll2002l . 

Historically, chemical and kinematic information were 
used as a basis to formulate the first galaxy formation mod- 
els. In their classical paper, lEggen. Lvnden-Bell fc Sandagd 
1 19621 ) analysed the properties and motion of 221 dwarfs 
and showed that those of lower metallicity tended to move 
on more highly eccentric orbits. The observed trends were 
interpreted as a signature that the stars now observed as 
a spheroidal halo formed during a rapid radial collapse 
that later co ntinued to form the stellar disk. About one 
decade later, ISearle fc ZinnI ([1973) measured the metallic- 
ities in a sample of globular clusters and found no signifi- 
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cant abundance gradient as a fu nction of the g alactocentric 
distance. These observations led lSearle fc Zinnl to formulate 
the hypothesis that the stellar halo formed over a longer 
timescale through agglomeration of many sub-galactic 'frag- 
ments' that may be similar to the surviving dwarf spheroidal 
galaxies today observed as satellites of the Milky Way. 
The observational evidence in support of this scenario has 
mounted significantly in the last decade ranging from detec- 
tion of significant dumpiness i n the phase space distribu- 



tion o f halo and disk stars (e.g. Maiewski. Munn fc Hawlev 
19961 : iHelmi et all Il999l : IChibafc BeersI 12000: Hcl mi et al 



20061 ) to the direct detection of satellite galaxies caught in 



ugnt in 
the act of tidal disruption (e.g. Ibata. Gilmore fc IrwirJ 19941 : 
iMartin et aTll2004l: IZucker et al.ll2006l). 



The ISearle fc ZinnI scenario, and the observational re- 
sults mentioned above, are in qualitative agreement with 
expectations from the hierarchical cold dark matter (CDM) 
model which, with the 'concordance' set of cosmological pa- 
rameters (ACDM), boasts a considerable degree of success 
in reproducing a large number of observations at various 
cosmic epochs, in particular on large scales. On galactic 
and sub-galactic scales, the success of the model has not 
been convincingly demonstrated yet, and a number of is- 
sues remain subject of a lively debate in the astronomical 
community. Among these are the dwarf galaxy counts, the 
degree of concentration of dark matter haloes, and the ro- 
tati on curves of dwarfs and low surface brightness galaxies 



(e.g. ISimon et al .1120031 : iKazantzidis et al.ll2004l : IStoehrll2006l : 
IStrigari et al.ll2007l . and references therein). 

The origin and structure of the stellar halo has 
been studied by several auth ors, u s ing a variety of 
techniques (for a review, see iHelmil |2008| ) . These in- 
clude cosmological numerical simulations with and with- 
out baryonic physics (e.g. iDi cmand, Madau & Moore 2003: 
lAbadi. Navarro fc Steinmeta 120061 ) , and phenomenological 
modelling of the evolution of baryons inside haloes, usually 
in combination with N-body si mulations that provide th e 
dynamical history of the system (jBullock fc Johnstonll2005r ) . 

A number of r ecent studies have raised concerns about 
the ISearle fc ZinnI scenario on the basis of the observation 
that stars in Local Group dwarf spheroidals (dSph's) tend 
to have lower a abundances than stars in the stellar hal o 
l|Shetrone et al.ll200ll : iTolstov et alll2003l : IVenn et al.ll2004 ) . 
The observed abundance pattern excludes the possibility 
that a significant contribution to the stellar halo comes from 
disrupted satellites similar to the present day dSph's. It is 
not entirely unexpected that the surviving satellites might 
be intrinsically different from the main contributors to the 
stellar halo, given they had a Hubble time to evolve as in- 
dependent entities. This argu ment has been put f orwar d 
in a recent serie s of papers bj^lBuUock fc Johnston! (|2005r ). 
[Robertson et al.l (|2005:;), and Font ct al. (2006) who repro- 
duce the observed chemical abundance pattern by combining 
mass accretion histories of galaxy-size haloes with a chemi- 
cal evolution model for individual satellites. In their model, 
the agreement with the observed trends is a consequence of 
the fact that the stellar halo originates from a few (rela- 
tively massive) satellites accreted early on and enriched in 
a-elements by type II supernovae. The surviving satellites 
are typically accreted much later, have more extended star 
formation histories and stellar populations enriched to solar 
level by both type II and type la supernovae. A more serious 



proble m with the Searle fc ZinnI scenario has been pointed 
out bv JHelmi et al.l ( 20061 ) who found a significant difference 
between the metal-poor tail of the dSph metallicity distribu- 
tion and that of the Galactic halo, demonstrating that the 
progenitors of present day dSph's are fundamentally differ- 
ent from the building blocks of our Galaxy, even at earliest 
epochs. One possible solution discussed in Helmi et al. is 
that the Galactic building blocks formed from the collapse 
of high-(T density fluctuations in the early Universe, while 
the pres ent day satellites would ori ginate from low-cr peaks 
(see also iSalvadori. Ferrara fc Schn eider 2008). 

The formation and the evolution of the baryonic com- 
ponent of galaxies is regulated by a number of non-linear 
processes operating on vastly different scales (e.g. shocking 
and cooling of gas, star formation, feedback by supernovae 
and active galactic nuclei, chemical enrichment, and stellar 
evolution). Most of these processes are quite poorly under- 
stood even when viewed in isolation. The difficulties grow 
considerably when one takes into account the fact that the 
physical properties of galaxies are determined by a com- 
plex network of actions, back-reactions, and self-regulation 
between the above mentioned physical processes. In recent 
years, different approaches have been developed to link the 
observed properties of luminous galaxies to the dark matter 
haloes within which they reside. Among these, semi-analytic 
models have developed into a powerful and widely used tool 
to study galaxy formation in the framework of the currently 
standard model for structure formation. In these models, the 
evolution of the baryonic component is modelled invoking 
'simple', yet physically and observationally motivated 'pre- 
scriptions'. These techniques find their seeds in the pioneer- 
ing work by W hite fc Reea (|l978l ): they have been laid out in 
a more detail ed form in the early 90s dW hitc fc Frcnk 199l|; 
[Colelll99ll : lKauffmann. White fc Guiderdonilll993 ) and have 
been substantially extended and refined in the last years by 
a number of different groups. Modern semi-analytic mod- 
els of galaxy formation take advantage of high resolution 
N-body simulations to specify the location and evolution of 
dark matter haloes - w hich are assumed to be the birthplaces 
of luminous galaxies (IKauffmann et al. 19991: [Benson et al.l 



I2OOOI : ISpringel et al.|[200ll : iHatton et al.ll2003i ). Using this 
'hybrid' approach, it is possible not only to predict observ- 
able physical properties such as luminosities, metallicities, 
star formation rates, etc., but also to provide full spatial 
and kinematical information of model galaxies thus allowing 
more accurate and straightforward comparisons with obser- 
vational data to be carried out. In this paper, we use this 
hybrid approach to study the formation of the Milky Way 
galaxy and of its stellar halo. 

The numerical simulations used in our study are de- 
scribed in Sec [2j while in Sec. |3] we give a brief introduc- 
tion to the adopted semi-analytic technique, and details of 
the specific model used in our study. In Sec. |4l we compare 
model results to the observed properties of the Milky Way, 
and discuss their dependence on a number of model parame- 
ters. In Sec.O we study the age and metallicity distribution 
of the spheroid and disk components of our model Galaxy. 
In Sec. [HI we analyse the formation and structure of the stel- 
lar halo. We discuss and summarise our results, and give our 
conclusions in Sec. [T] 
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Table 1. Numerical parameters of the re-simulations used in this 
study. In the table, we give the particle mass rrip, the number of 
particles, the starting redshift ^atart, and the gravitational soft- 
ening e. All these quantities refer to the particles in the high 
resolution region only. 



Table 2. M200 and Vi, 
study. 



for the four simulations used in this 



Name nip [fe'^M©] 



Nv 



e [h-^kpc] 



GAO 1.677 X 10* 

GAl 1.796 X 10^ 

GA2 1.925 X 10^ 

GAS 2.063 X 10^ 



68323 


70 


1.4 


637966 


80 


0.8 


5953033 


90 


0.38 


55564205 


60 


0.18 



2 THE SIMULATIONS 

In this paper, we use the re-simu lations of a 'Milky Way' 
halo (the GA serie s) described in IStoehr et al.l (|2002l ) and 
IStoehr et al.l ((2003) • The underlying cosmological model is 
a flat A-dominated CDM Universe with cosmological pa- 
rameters: fim = 0.3, f^A = 0.7, Ho = 70kms~^Mpc"\ 
n = 1, and as — 0.9. Th e simulations were generated us - 
ing the 'zoom' technique (|Tormen. Bouchet fc Whitelll997l ) 
starting from an intermediate-resolution simulation (parti- 
cle mass ~ 10* Mq) of a 'typical' region of the Universe. A 
'Milky Way' halo was selected as a relatively isolated halo 
which suffered its last major merger at 2: > 2, and with ap- 
proximately the correct peak rotation velocity. These choices 
were intended to select a candidate for resimulations with 
properties that closely match the observed ones: our Galaxy 
is believed to have had a quiet merging history, with the 
last major merging event occurring some 10-12 Gyr ago 
l|Gilmore et al.ll2002r ). and it is a quite isolated galaxy with 
the closest cluster (Virgo) lying at a distance of about 20 
Mpc. 

It is important to note that haloes with such quiet merg- 
ing histories are not uncommo n in ACDM cosmolog ies. Us- 
ing the Millennium Simulation (ISpringel et al.ll2005l ) , we de- 
termined the epoch of the last accretion event correspond- 
ing to an increase in mass of the main progenitor of a factor 
~ 0.3 (i.e. a 'major merger') for haloes with present-day 
mass in the range 1 to 3 x lO^'^ Mq . We find that about 57 
per cent of these haloes experienced their last 'major merger' 
at redshift larger than 2. The estimated fraction should be 
considered as an upper limit, as part of the accretion might 
occur through diffuse material. Lower - but still significant 
- fractions ( ^ about 45 per c ent) h ave been found in a recent 
analysis bv lStewart et al.l (|2007n . and in other independent 
studies (Boylan-Kolchin, private communication). 

The selected halo was then re-simulated at a series of 
four progressive ly higher resolution simulations using the 
code GADGET ijSpringel. Yoshida fc Whitell200lh . The nu- 
merical parameters used in the high resolution regions of 
these re-simulations are summarised in Table [T] 

Simulation data were stored in 108 outputs from z — 
37.6 to z = 0. These are approximately logarithmically 
spaced in time down to z — 2.3, and approximately linearly 
spaced in time thereafter. For each simulation snapshot, we 
constructed group catalogues using a standard friends-of- 
friends (FOF) algorithm with a linking length of 0.2 in units 
of the mean particle separation. Each group was then de- 
composed into a set of disjoint subst ructures using the algo- 
rithm SUBFIND (jSpringel et alllioOli ) . As in previous work, 



Name M200 [Mq] Vmax [kms" 

GAO 3.35 X 10^2 250.71 

GAl 3.28 X 10^2 246.60 

GA2 3.21 X 10^2 247.05 

GA3 2.98 X 10^2 251.20 



we consider to be genuine substructures only those with at 
least 20 bound particles which sets the subhalo detection 



l9 



,7 



de- 



limit to 4.79 X Iff', 5.13 X 10", 5.5 x 10', and 5.89 x 10'' Mq 
for the four simulations used in our study. Substructure cat- 
alogues were then used to construct merger h istory trees for 
all self-bou nd haloes as described in de tail in lSpringel et al.l 
(2005) and iDe Lucia fc Blaizo^ (|2007l ). We recall that the 
merger tree construction is based on the determination of 
a unique descendant for any given halo. In order to deter- 
mine the appropriate descendant, for each halo we find all 
haloes in the following snapshot that contain its particles, 
and then count the particles giving higher weight to those 
that are more t ightly bound in the h alo under consideration. 
As a reference, ISpringel et al.l l|2005l ) gave about three times 
more weight to the half most bound particles of a halo of 
100 particles. (The relative weight of the half most bound 
particles increases with particle number.) For our analysis, 
we are interested in tracing well the cores of the accreted 
satellites (see Sec. [5}, and we noted that the tree construc- 
tion adopted in previous work lead to occasional premature 
mergers. In order to avoid these events, we have slightly in- 
creased the weight of the most bound particles (four times 
more weight is given to the half most bound particles of a 
100 particle halo). The merger trees constructed as described 
above, represent the basic input needed for the semi-analytic 
model which is described in the next section. 

In Table [2j we list the masses (M200) and maxi- 
mum velocities for the main halo in the four simula- 
tions used in this study, [j Latest observational results give 
Vmax ^ 180 - 220km s ~^ and a mass Mm w ~ 1 x 10^^ M© 
(JBattaglia et all l2005l : ISmith et all |2007| ). Our simulated 
haloes therefor e appear mo r e mas sive than the Milky- Way 
halo. Following iHelmi et al.l l|2003r ). we scale our simulations 
to a 'Milky Way' halo by adopting a scaling factor in mass 
M200/MMW = 7^ = 2.86. This implies that we scale down 
the positions and velocities by the same factor 7 = 1.42 in 
all four simulations. 

Fig. \T\ shows the projected density distribution of the 
dark matter for the simulation GA3, at six different red- 
shifts. The box at 2; = is centred on the most bound par- 
ticle of the GA3 halo, while those at higher redshift are 
centred on the most bound particle of the main progenitor 
of the GA3 halo at the corresponding redshift. 



M200 is defined here as the mass within a sphere of density 200 
times the critical density. 
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Figure 1. Projected density distribution of the dark matter for the simulation GAS, at six different redshifts (from top left to bottom 
right: 10.4, 4.8,2.0,1.0,0.5, and 0). The box at z = is centred on the most bound particle of the GAS halo, while those at higher redshift 
are centred on the most bound particle of the main progenitor of the GAS halo at the corresponding redshift. For all panels, the depth 
of the box is 500 kpc comoving. As explained in the text, positions are scaled by a factor 7 = 1.42. 
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3 THE HYBRID MODEL OF GALAXY 
FORMATION 

A complete review of semi-analytic techniques goes beyond 
the scope of this paper. For a pedagogical and historical in- 
troduction to semi-analytic m ethods, we ref er the interested 
reader to the recent review bv lBaughl l|2006r ). In this section, 
we briefly summarise how a semi-analytic model is grafted 
onto high resolution A^-body simulations, and give a brief 
account of those aspects of the physical model which are 
relevant for the present study. 

The basic assumption is that galaxies form when gas 
condenses at the centre of dark matter haloes. Star forma- 
tion, feedback processes, chemical enrichment, etc. take then 
place according to analytical laws which are based on the- 
oretical and/or observational arguments. Adopting this for- 
malism, it is possible to express the full process of galaxy 
evolution through a set of differential equations that describe 
the variation in mass as a function of time of the different 
galactic components (e.g. stars, gas, metals), and that are 
coupled to the merger history of the dark matter haloes ex- 
tracted fr om the A'^-bo dy simulations (e.g. see Fig. 1 and 
Sec. 4.7 in lDe Lucia. K auffmann & White 2004). Given our 
limited understanding of the physical processes that regu- 
late galaxy formation and evolution, the equations describ- 
ing these processes contain 'free' parameters whose value is 
typically chosen in order to provide a reasonably good agree- 
ment with the observational data in the local Universe. 

The semi-analytic model used in this s tudy builds on 
the m etho dology originally intro duced by ISpringel et al.l 
l|200ll ) and lOe Lucia et all (|2004l '). and has been recently 
updated to include a model for the suppression of cooling 
fl ows by 'radio-mode' AGN feedback as described in detail 
in lCroton et al.l (|2006l ). Details about the modelling adopted 
fo r the various physic al pr ocesses considered can be found 
in lCroton eTaH l|2006l ) and lDe Lucia fc Blaizo3 l|2007h . 

We recall that our approach follows explicitly dark mat- 
ter haloes when they are accreted into larger systems. This 
allows the dynamics of satellite galaxies residing in the in- 
falling haloes to be properly followed until their parent dark 
matter substructures a re completely destroyed by tidal trun- 
catio n and stripping IjGhigna et al.l l2000l : iDe Lucia et al.l 
l2004l : iGao et aLrl2004h . When this happens, the satellite 
galaxy residing at the centre of the substructure under con- 
sideration is assumed to merge onto the central galaxy after 
a residual surviving time that is estimated from the rela- 
tive orbit of the two merging objects, at the time of subhalo 
disruption: 



1.17 ■/< 



Kir 



fudge 



In A ■ G • Ms, 



In the above equation, Mgat is the mass of the substruc- 
ture at the last time it is identified, rsat is the distance 
between the merging halo and the centre of the structure 
on which it is accreted, K/iriai is the virial velocity of the 
accreting structure, and we approximate the Coulomb loga- 
rithm with InA = (1 -|- Mvir/Afaat). We therefore allow the 
satellite galaxy to further sink in through dynamical fric- 
tion, even after its dark matter halo has fallen below the 
resolution limit of our simulations. When a halo is accreted 
onto a larger structure, its mass is reduced by tidal strip- 
ping and its orbit is shrunk by dynamical friction, until the 



substructure can no longer be identified as a self-bound over- 
density orbiting the smooth dark matter background of the 
larger system. This 'disruption' of the dark matter substruc- 
ture typically occurs at rsat ^ 0.1-Rvir which can be much 
larger than the separation from which the galaxy merger 
is expected to happen. We note, finally, that previous work 
has shown that this treatment of orphan galaxies is required 
in order to reproduce the correlation signal a t small scales 
(|Wang et al.ll2006l : lKitzbichler fc Whitell2008D . 

Mild variations of the above formula are used in dif- 
ferent semi-analytic models with differences entering mainly 
the way each model treats the orbital distribution. In A, and 
■f fudge (see discussion in iBovlan-Kolchin. Ma fc Quataert 
l2008h . In this work, we have followed iDe Lucia fc Blaizot 
(2003) and used /fudge = 2, which is in better agree- 
ment with recent numerical work indicating that the 
classical dynamical friction formulation tends to under- 
estimate the merging times measured from simulations 
l|Bovlan-Kolchin et al.l boOSi : I Jiang et aLlbOOSl ). Results of 
these studies, however, are still quantitatively different and 
additional work is needed in order to calibrate the correct 
pre-factor and/or corrections to apply to the classical for- 
mula. 

The stellar mass of the satellite galaxies is assumed to 
be unaffected by the tidal stripping process that reduces 
the mass of its parent halo. However, such tidal effects may 
well dominate the disruption process of sate llites orbiting in 
a galaxy halo IjPefiarrubia fc Bensonll2005h . Therefore, our 
implementation could lead to an underestimate of the stellar 
mass in diffuse form, and to an overestimate of the number 
and/or luminosity of satellite galaxies. We will see however, 
that this does not affect our conclusions in any significant 

way. 

As in previous work (jCroton et al.ll2006l : iDe Lucia et al.l 

120061 ) ■ we assume that spheroid formation occurs through 
both mergers and disk instability. In the case of a 'minor' 
merger, we transfer the stellar mass of the merged galaxy 
to the spheroid of the central galaxy. The cold gas of the 
satellite galaxy is added to the disk of the central galaxy, 
and a fraction of the combined cold gas from both galax- 
ies is turned into stars as a result of the merger. Any stars 
that formed during the burst are also added to the disk of 
the central galaxy. The photometric properties of the galaxy 
ar e updated accordingly using the same method described 
in iDe Lucia et al.l (|2004l ). If the mass ratio of the merging 
galaxies is larger than 0.3, we assume that we witness a 'ma- 
jor' merger that gives rise to a more significant starburst and 
destroys the disk of the central galaxy completely, leaving 
a purely spheroidal stellar remnant. The galaxy can grow 
a new disk later on, provided it is fed by an appreciable 
cooling fiow. Note that the modelled spheroidal component 
includes both the bulge and the s t ellar h alo. 

Following [Mo. Mao fc White! (|l998h . we assume that a 
stellar disk becomes unstable when the following condition 
is verified: 



Vn. 



(Gmdisk/'-diak)^''" 



< 1.1 



This condition is based on the numerical work by 
lEfstathiou. Lake fc Negropontd l) 19821 ') who used A''-body 
simulations to investigate the development of global insta- 
bilities in exponential disks embedded in a variety of haloes. 
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In this work, we assume that each time a galaxy meets the 
above criterion, a fixed fraction (-Finat) of the stellar mass 
in the disk is transferred to the central spheroid compo- 
nent. The instability criterion is applied only to disk dom- 
inated systems (Msphcroid/Mtot < 0.1). This choice is mo- 
tivated by numerical studies that showed that a spheroid 
component can stabilise the disk against bar formation, 
giving the galaxy an inner Lindblad resonance which does 
not allow swing amplification of waves thro ugh the centre 
fe.g. lSellwoo"dlll989l : [Seriwood fc Moordl 19991 ). The parame- 
ter -Finst is chosen so as to get a morphological mix that is in 
reasonable agreement with the observational measurements 
for the local Universqj, and is assumed to be equal to 0.5. 
We note that, according to the condition above, the occur- 
rence of disk instability episodes follows essentially from the 
way the disk mass builds up relative to the total mass of the 
halo, and that the condition Msphcroid/Mtot < 0.1 efficiently 
suppresses the occurrence of late instability episodes even in 
galaxies with modest bulge-to-total ratios. 

Admittedly, our modelling of disk instability is very 
simplified, and gives a very crude description of the 
complex phenomenology associated to bar formation and 
evolutio n, which is the subject of c ur rent active re- 
search (IMartinez-Valpuesta et al.l 120061 : iBerentzen et al.l 
120071 : ICurir et al.ll2007l . and references therein). We are, for 
example, neglecting the possibility that bar formation pro- 
duces an inflow of gas towards the centre that could fuel 
starburst/AGN activity, and that can eventually lead to bar 
disruption. In addition, present simulations do not provide 
clear indications about the fraction of disk mass that gets re- 
distributed and how this depends on the halo/galaxy prop- 
erties. We therefore consider the results of our simple mod- 
elling as just indicative and, in the following, we also com- 
ment on results obtained when the disk instability channel 
for bulge formation is switched off. 

In this paper, we adopt an improved model to estimate 
the disk radii. Assuming conservation of specific angular mo- 
mentum, we assume that when hot gas cools at the centre 
of dark matter haloes, it settles in a rotationally supported 
disk with exponential scale-length given by: 

where A is the halo spin p arameter IjMo et aLlllQQa ). Fol- 
lowing |Hatton_eLalJ ((2003), we recompute the scale-length 
at each time-step by taking the mass-weighted average gas 
profile of the disk and that of the new material being ac- 
creted. The scale-length of the disk is not altered after a 
galaxy becomes a satellite. 

As in previous work, we assume that the star formation 
occurs at a rate given by: 

ifj = aSpMsf/tdyn 

where idyn = ''disk/Vvir is the dynamical time of the galaxy 
and we assume that the star forming region (rdisk) extends 
to about 3 X ro. The parameter asp regulates the efficiency 
of the conversion of gas into stars, and M^t represents the 
amount of gas available for star formation which can be 
expressed as: 



^ The morphological mix resulting from the adopted model has 
been computed for the Millennium Simulation. 



Ms) 



2-K 



T,D{r)rdr 



where rcrit is the radius at which the gas surface dens ity 
drops below the following critical value l|Kennicuttlll989l ): 

Scrit[M0pc-2] = 0.591/[kms-V^disk[kpc]. 

For a complete summary of the model param eters adopted 
in our fiducial model, we refer to Table 1 of 
(2009). As explained in iDe Lucia fc Blaizot 



Croton et al.l 
(20071), 



the 

adoption of a Chabrier Initial Mass Function led to slight 
modifications of some of these parameters. We refer to these 
papers for more details about the parametrisation adopted 
for the various physical processes explicitly taken into ac- 
count in our model, and for a complete characterisation of 
model parameters. 



4 DEPENDENCY ON MODEL PARAMETERS 
AND NUMERICAL RESOLUTION 

In order to provide an illustration of how model results are 
affected by different choices of model parameters, we show in 
Fig. [2] how different physical properties of our model Milky 
Way galaxy vary on a limited model grid where we have al- 
tered only the values of the star formation and supernovae 
feedback efficiencies. These corr espond to the param eters 
"asp" and "edisk" in Table 1 of ICroton et all (|2006l) . and 
regulate the amount of cold gas that is converted into stars 
in a disk dynamical time, and the amount of cold gas re- 
heated by supernovae explosions. As explained in the previ- 
ous section, our fiducial model uses the same combinatio n 
of model parameters adopted in iDe Lucia fc BlaizotI ((20071) , 
and is indicated by the vertical solid line in Fig. [21 

In the figure, different colours for filled circles corre- 
spond to different values for the star formation efficiency 
parameter: 0.03 (red), 0.08 (green), 0.13 (blue). Different 
colours for open circles correspond to different values for 
the SN feedback efficiency: 1.5 (green), 3.5 (red), 5.5 (blue). 
Red horizontal lines in each panel indicate the observational 
measurements. Fig. [2] shows results from our highest resolu- 
tion simulation (GA3), but similar trends are obtained for 
the lower resolution simulations used in this study (see be- 
low). 

Increasing qsp (red, green, blue filled circles) and keep- 
ing Edisk constant (open circles of the same colour) pro- 
duces an increase of total mass, and a corresponding de- 
cr ease in of the amount of cold gas available. As explained 
in lDe Lucia et al.l (I2OOJ), metals are exchanged between dif- 
ferent components proportionally to the exchanged mass so 
that the increase in mass is reffected in a parallel increase 
in the amount of metals (bottom left panel in Fig. [21). The 
mass of the spheroid component (top left panel) varies in a 
non-monotonic way with increasing star formation efficiency 
because disk instability episodes (see previous section) occur 
at different times. In our model, black holes grow primar- 
ily during mergers, both by merging with each other and 
by accretion of cold gas (see Sec. 3.4 of Croton ct al. 2006)). 
The black hole mass therefore decreases with increasing star 
formation efficiency (middle left panel) as a consequence of 
the decrease of the cold gas available. 

Increasing edisk (green, red, blue open circles) and keep- 
ing asp constant (filled circles of the same colour), the final 
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Figure 2. Physical properties of our model Milky Way galaxy from the simulation GAS, as a function of different combinations of the 
star formation efficiency (q:sf)i a-nd of the supernovae feedback efficiency (edisk)- Different colours for filled circles correspond to different 
values for the star formation efficiency parameter: 0.03 (red), 0.08 (green), 0.13 (blue). Different colours for open circles correspond to 
different values for the SN feedback efficiency; 1.5 (green), 3.5 (red), 5.5 (blue). The vertical black line (corresponding to model number 
4) indicates our fiducial model. Red horizontal lines in each panel indicate observational estimates. Note that in the bottom right panel 
we compare the metallicity of the spheroid in our models to that of the Galactic bulge. 



amount of gas stays almost constant because gas ejected by 
satellite galaxies is rapidly re-incorporateqj, and because the 
higher level of enrichment of the hot gas component causes 
a more efficient cooling. As a consequence, the total stel- 
lar mass increases (top left panel), and the black hole mass 
decreases slightly (middle left panel). 

For our reference model, the total stellar mass is ~ 
6 X 10^° Mq, in very good agreement with the estimated 
value of ~ 5 — 8 X 10^" Mq. The mass of the spheroidal 
component is instead slightly lower than the observed value 
(assumed to be ~ 25 p er cent of the disk s tellar mass, or 



1.25 - 1.6 X 10^° Mq - iBissantz et al.ll2004l). Studymg the 



proper motion of Galactic Centre stars. ISchodel et ahl (|2002l ) 
estimated a central point mass of (3.7 ± 1.5) x 10^ Mq, 
which is lower than the value obtained for our Milky Way 



^ Following iDe Lucia et al.l ll2004h . we assume that the material 
that is ejected outside the haloes as a consequence of supernovae 
explosions, can be re-incorporated into the hot gas component 
after a time-scale that is proportional to the dynamical time-scale 
of the halo. O ur default model ado pts the same re-incorporation 
efficiency as in lCroton et al.l 1120061 ). 



in the reference model. This is not surprising given that our 
model for the formation of black holes was tuned to match 
the local relation b etween black hole mass and bulge mass 
IjCroton et al.l 120061 ). and that the value measured for the 
Milky Way is offset from this relation. The total mass of 
H2 in the Milky Way galaxy is ~ 1 x 10^ M^^^ and the total 
mass of HI gas is ~ 5 x 10^ Mq (iBlitjIlQQTT and references 
therein). Fig. [2] shows that our fiducial model gives a gas 
mass which is about twice the estimated value, and that this 
'problem' could be alleviated by increasing the star forma- 
tion efficiency. The stars of the Galactic bulge are found to 
peak at near solar value ([Fe/H]~ —0.15), with a relatively 
wide spread, while the stars in the disk have [Fe/H]~ —0.1 
IjFreeman fc Bland-Hawthorri2002l . see also Sec. [5}. The cor- 
responding metallicities for the galaxy in the reference model 
are not far from the observed values. Note, however, that the 
metallicities shown in Fig. |2] and in the rest of this paper, 
refer to total metallicities (and not only to [Fe/H]). It should 
be stressed that our model assumes an instantaneous recy- 
cling approximation, i.e. we neglect the delay between star 
formation and the recycling of gas and metals from stellar 
winds and supernovae. The model is therefore not able to 
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take into account the evolution of different element abun- 
dances, and in particular it does not describe well the ele- 
ments around the iron-peak, which are mainly produced by 
supernovae la and which typically return their products on 
longer timescales. We will come back to this issue later. 

Overall, Fig. [2] demonstrates that our reference model 
is in relatively good agreement with observational measure- 
ments (although there are other combinations of model pa- 
rameters that also provide a reasonable agreement). This 
is not entirely surprising given that model parameters are 
tuned to reproduce the local galaxy luminosity function 
and the mass and lu minosity of 'Milky- Way' galaxies (see 
IPe Lucia et al.l 12004 1 . In previous work, we have shown 
that our fiducial model is in quite good agreement with 
the observe d relations b etween stellar mass, gas mass, and 
metallicity iJDe Lucia et al. 2004), the observed luminosity , 
colour, and morphol ogy distributions (ICroton et al.l 120061 : 
iDe Lucia et al.l 120061) . and the observed two-point corre - 
lation functions llSpringel et al.1 l2005l : IWang et al] l2007l l. 
iKitzbichler fc White! ([20071) have recently shown that it also 
agrees reasonably well with the observed galaxy luminosity 
and mass functions at higher redshift. In the following, we 
will therefore focus on results from our fiducial model. It 
is worth mentioning at this point that our d efault model, 
which assumes a redshift of reionization z ~ 7 (jCroton et al.l 
120061 ). over-predicts the number of observed satellites. The 
situation is significantly improved when adopting an earlier 
epoch of reionization which does not affect significantly the 
model Galaxy and the results presented in this paper. We 
will present a more detailed analysis of the satellite popula- 
tion for our model Milky Way in a forthcoming paper. 

Fig. [3] shows the evolution of different mass and metal- 
licity components for the model Milky Way galaxies in 
the four simulations used in this study (lines of different 
colours). The histories shown in panels (a) to (f) are ob- 
tained by tracking the evolution of the 'main progenitor', 
which is obtained by linking the galaxy at each time-step to 
the progenitor with the largest stellar mass. For our model 
Milky Way, the galaxies merging onto the main branch have 
stellar masses that are much smaller than the current mass 
of the main progenitor, over most of the galaxy's life-time, 
explaining the smooth increase of the stellar mass compo- 
nent. In these cases, the main progenitor branch does indeed 
provide a quite good approxi mation of the evolution of the 
galaxy itself (see discussion in lDe Lucia fc Blaizotll2007^ . 

Fig. |3] shows that approximately half of the final mass 
in the dark matter halo is already in place (in the main 
progenitor) at 2: ~ 1.2 (panel a) while about half of the 
final total (disk + spheroid) stellar mass is only in place at 
z ~ 0.8 (panel b). About 20 per cent of this stellar mass 
is already in a spheroidal component (panel c). The mass 
of the spheroidal component grows in discrete steps as a 
consequence of our assumption that it grows during mergers 
and disk instability events, and approximately half of its 
final mass is already in a spheroidal component at z ~ 2.5. 
The cold gas mass contenl|J, in contrast, varies much more 
gradually (panel d). Fig. |3] also shows that there is a very 



^ Note that in our model, there is no cold gas in the spheroidal 
component. So all gas in panel (d) of Fig. |3]is associated to the 
disk. 
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Figure 4. Star formation history of the model Milky Way galax- 
ies from the four simulations used in this study (different colours) . 



weak increase of the stellar metallicity at late times: it only 
varies from ~ 0.75 solar to ~ 0.9 solar over the last 5 Gyr 
(panel e), while the gas-phase metallicity varies from ~ 1 
solar to ~ 1.3 solar over the same interval of time. The 
more efficient enrichment of the gas component is due to 
the fact that in our model all metals produced by new stars 
are instantaneously returned to the cold phase, i.e. we are 
assuming a 100 per cent mixing efficiency of the metals with 
the cold gas already present at the time of star formation. 

Interestingly, our model produces consistent evolutions 
for all four simulations used in our study, despite a very 
large increase in numerical resolution (see Table [l]). This 
is particularly true for the stellar mass (panel a) and for 
the stellar metallicity (panel e). The amounts of gas (panel 
d) and metals in the gas phase (panel f) exhibit a much 
more noisy behaviour, but the overall evolution is still very 
similar. We note, however, that some panels (e.g. panels b, 
e, and f) do not show convergence of the results and this is 
driven by the lack of complete convergence in the A'^-body 
simulations (panel a shows a clear difference between GAS 
and the lowest resolution simulations). 

Fig. 13] shows the star formation history for the model 
Milky Way galaxies from the four simulations used in this 
study (colour coding as in Fig. [S}. The star formation his- 
tories shown in Fig. 3] exhibit a quite bursty behaviour 
with several broad periods of enhanced star formation. In 
our models, these bursts are triggered by minor merger 
events. This intermittent behaviour is in qualitative agree- 
ment with the trend recovere d by chromospheric age dis- 
tributions of late-type dwarfs (|Rocha-Pinto et al.ll2000l . and 
references therein). Fig. 2] also shows that the star forma- 
tion rate declined slightly over the last ~ 4 Gyr and has a 
current value of ~ 5MQyr~^. This trend is also in qual- 
itative agreement with observational measurements (e.g. 
Hernandez. Valls-Gabaud fc Gilmor^ I2OO0I : iBertelli fc NeSl 
200l|), but the level of star formation activity is significantly 
higher than the < IMQyr"^ estimated for the mean star 
formation rate in the Milky Way disc over the last few Gyrs 
(|Just fc Jahreis3i2007l ). Although the general behaviour is 
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very similar for all four simulations used in this study, there 
is again no perfect convergence of the results, reflecting the 
behaviour of the available cold gas (panel d in Fig. [Sj . 

For our default model and for the four simulations used 
in this study, the spheroid mass varies between 6.4 x 10^ M© 
and 7.3 x 10^ M© which is slightly lower than but still com- 
patible with the observational estimate. As explained in the 
previous section, our model includes different channels for 
the formation of a spheroidal component: mergers (both mi- 
nor and major) and disk instability. Fig. [3] (panel c) shows 
how the total mass in the spheroidal component grows as 



a function of redshift. In Fig. [5] we show the cumulative 
fraction of the mass in the spheroidal component coming 
from minor mergers (solid coloured lines), and disk insta- 
bility (dashed lines). Our model Milky Way galaxies do not 
experience any major merger during their life-times. In all 
four simulations used in our study, there are two episodes 
of disk instability, the flrst ~ 12 Gyr ago and the second 
~ 1 Gyr later. Disk instability contributes in the range of 
48 (GAl) to 58 (GA2) per cent of the final steUar mass 
in the spheroidal component. Minor mergers contribute to 
the remaining stellar mass in the spheroidal component, and 



© 2008 RAS, MNRAS OOO.ITH201 



10 G. De Lucia & A. Helmi 



0.8 
0.6 


GAO . 
GAl . 
GA2 _ 




, ^^g 










. 4 

















0.4 




-| 


1 


^ 


0.2 


- 




1 







L__j=v — \ ' J 


on 


,^^^^.^1:-, , 



4 6 8 10 

lookbacktime (Gyr) 



12 14 



Figure 5. Cumulative fraction of tlie mass in the spheroidal com- 
ponent coming from minor mergers (solid lines) and from disk 
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all occur at later times with respect to the disk instability 
episodes. As explained in Sec. [S] the stars formed during 
the bursts accompanying minor merger events, are added to 
the disk component. As a consequence, no spheroid star is 
formed in situ in our models. 

As noted in Sec.|3l our model for disk instability is quite 
simplified and, given the uncertainties involved, the fraction 
of mass contributed by this channel should be considered 
as indicative. The results shown in Fig. [5] and Fig. [2] indi- 
cate that our default model would not produce a spheroid 
that is as massive as observed if the disk instability channel 
is switched off. On the other hand, it is now fairl y well es- 
tablis hed that the Milky Way is a barred galaxy l) Gerhard! 
12002 . and references therein). Since global disk instabilities 
are commonly considered the main mechanism for the for- 
mation of a bar, the presence of a bar in the Milky Way is 
a direct indication that it has experienced episodes of disk 
instability during its life-time. 



5 METALLICITY AND AGE DISTRIBUTIONS 

As mentioned in Sec. [l] accurate measurements of ages, 
metallicity and kinematics have been collected for a large 
number of individual stars of our Galaxy, and a large amount 
of new observational measurements are expected in the near 
future. The available information is n icely summarised in 
Fig. 2 of iFreeman fc Bland-HawthornI l|2002l ). 

The available data are indicative of a significant pop- 
ulation of stars in the local disk with ages < 8 Gyr, 
and metallicity distribution (mostly from G/K dwarfs) 
that peaks at [Fe/H]~ —0.1. From an analysis of the 
main sequence turn-off stars in the Hipparcos dataset, 
iBinnev. Dehnen fc Bertellil (120001 ) obtained a best fit age for 
the oldest disk stars of > 11 Gyr. The properties of the thick 
disk are relatively well known only within a few kpc from the 
Sun. In this region, typical thick disk stars have ages as old 



as ~ 12 Gyr and intermediate metallicities ([Fc/H]~ —0.6). 
For the Galactic bulge, deep HST and ISO colour-magnitude 
diagrams suggest a dominant old (> 10 Gyr) population 
(JFeltzing fc Gilmordlioool : Ivan Loon et "all |2003| ). The ISO 
data also suggest the presence of a small intermediate-ag e 
component which is traced by OH/IR stars l|Sevensterll999l ). 
but the interpretation is complicated by thin disk contami- 
nation. Bulge stars have a metallicity distribution that peaks 
at [Fe/H]~ — 0.15 dex, with a broad range and a tail to low 
abundances dM cWilliam fc Rich' Il994 IZoccali et all l200l 
iFulbright. McW illiam & Rich 200^). 

Fig. [6] shows the age and metallicity distribution of all 
stars (left panel) and of the stars in the spheroidal com- 
ponent (right panel) of our model Milky Way galaxy from 
the highest resolution simulation used in this study (results 
for lower resolution simulations are very similar). In each 
panel, black and white contours enclose the regions con- 
taining 99.6, and 60 per cent of the stars. On the right 
panel, solid contours refer to our fiducial model where stars 
can be transferred from the stellar disk component to the 
spheroid through disk instability. Dashed lines refer to a 
model with the same combination of parameters but with 
the disk instability channel switched off. Fig. [S] shows that 
the Galaxy contains stars of all ages - which reflects the pro- 
longed star formation history shown in panel (g) of Fig. [Bl- 
and covering a relatively limited range of metallicities cen- 
tred around solar. The stars in the spheroidal component 
all have very old ages (> 11 Gyr), their metallicity distri- 
bution peaks at sub-solar value and exhibit a pronounced 
tail of low metallicity stars. If spheroid growth through 
disk instability is switched off, the stars in the spheroidal 
component have a broader tail towards low and high abun- 
dances. The left panel of Fig. [6] shows a very shallow age- 
metallicity relation with a broader tail towards lower metal- 
licities for older stars, in qualitative agreement with ob- 
servational measurements (INordstrom et al.l |2004| . see also 
iHolmberg. Nordstrom fc AndersenI 120071 ) . We caution the 
reader that the results shown in Fig. [S] represent a global 
average while observed samples represent a 'local sample'. 
In addition, the model age-metallicity distributions are not 
convolved with typical observational errors which naturally 
tend to broaden the distributions. 

The metallicity distributions (i.e. the projections along 
the y-axis of the maps shown in Fig. [6l of the stars in 
the disk and spheroid of our model galaxies are shown in 
Fig. [T] The left panels show the metallicity distribution of 
all stars (black histograms) and of the stars in the disk 
(blue histograms) compared to t he observational measure- 
ments bv iWvse fc Gilmord (119951) (ora nge histogram in the 
top panel) and INordstrom et al.l (12004) (green histog r am in 
the bottom panel). The sample of IWvse fc Gilmord (|l995r ) 
is a vo lume sample of long-liv ed thin disk G stars, while 
that of lNordstrom et all (|2004l ) contains a large (~ 14, 000) 
number of F/G dwarfs in the Solar neighbourhood (located 
at distances < 100 pc). In both these studies, metallici- 
ties were derived using a photometric calibration based on 
a relationship between [Fe/H] and Stromgren colours (e.g. 
ISchuster fc Nissenlll989l ). 

Fig. [7] thus shows that the metallicity distribution of 
disk stars in our model Milky Way peaks at approximately 
the same value as observed, but that it exhibits a defi- 
ciency of low metallicity stars. In comparing the model 
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and observed metallicity distributions, however, two factors 
should be taken into account: (1) the observational metal- 
licity measurements have some uncertainties (~ 0.2 dex) 
which tend to broaden the true underlying distribution; (2) 
the metallicity estimate used is an indicator of the iron 
abundance, which is mainly produced by supernovae la 
and therefore not well described by the instantaneous re- 
cycling approximation adopted in this study. In order to 
estimate the importance of this second caveat we have con- 
verted the measured [Fe/H] into [0/H] using a linear rela- 
tion, which is obtained by fitting data for thin disk stars 
from iBensbv. Feltzing fc LundstrSmI l|2004l ). The result of 
this conversion is shown by the dashed orange histogram 
in the top panel of Fig. [7] for the original measurements by 
IWvse fc Gilmord . The observed [0/H] metallicity distribu- 
tion is now much closer to the modelled log[Z/Z0] distribu- 
tion. 

The right panels in Fig. [7] show the metallicity distri- 
bution of the spheroid stars in our fiducial model (red his- 
togram) and in a model where the disk instability channel is 
switched off (dashed black histograms). In what follows, it 
may be useful to associate the stars originating in the disk 
instability to the bulge, and the rest of the spheroid stars to 
the stellar halo. Model results a r e com pared to observational 
measurements bv lZoccali et al.l l|2003l) (orange histogram in 
the top panel) an d recen t results bv iLecureur et al.l (|2007l ) 
and IZoccali et al.l ([2003) based on a sample of about 400 
K-giants in two Galactic bulge windows, all observed with 
high dispersion spectra (R > 20000) with GIRAFFE on 
VLT (green histogram in the bottom panel). These panels 
show that the metallicity distribution of the model spheroid 
peaks at lower value than observed, and that it contains a 



smaller fraction of high metallicity stars. Comparison be- 
tween the solid red histogram and the black dashed his- 
togram shows that the disk instability is responsible for 
the pronounced peak around Log[Z/ZQ]~ —0.25 due to the 
transfer of a large fraction of disk stars into the spheroidal 
component, and which may be associated to the Galactic 
bulge component (Fig. [3] and Fig. [5] show that the mean 
metallicity of the stellar component is ~ 0.6 Zq at the time 
of the major episode of disk instability). The same caveat 
discussed above applies to the metallicity distribution of 
the spheroidal components. We note that a conversion from 
[Fe/H] to [O/H] of the observed metallicity scale would bring 
most of the observed bulge stars to [0/H]> —0.2 suggest- 
ing that our model spheroid is significantly less metal rich 
than the observ ed Ga lactic bulge. However, we note that 
Fulbright et al.| l|2006l ) find a slightly more metal-poor dis- 
tribution for the 'Baade's Window', c ompared to t h e dis- 
tribution found for the same field by IZoccali et al.l (|2008l . 
plotted in the bottom right panel of Fig. [7]). I n addition, an 
outer bulge field (located at fe = —12 deg) from lZoccali et al.l 
( 2008 ) peaks at lower metallicity compared to the metallic- 
ity distribution from the other two fields presented in the 
same study. This suggests that the uncertainties (as well as 
field to field variations) in the metallicity distribution of the 
bulge may alleviate the disagreement with our model results. 



The metallicity distribution of the remainder spheroid 
stars (i.e. those originating in minor mergers) is relatively 
broad, and more metal-rich than observed for the Galactic 
stellar halo. 
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Figure 7. Motallicity distribution for stars in the disk (blue histograms in the left panels) and spheroid (red histograms in the right 
panels) of the model Milky Way from the highest resolution simulation in this study (GAS). The solid black histograms in the left 
panels show the metallicity distribution for all stars in the model galaxy, while the dashed black histograms in the right panels show 
the metallicity distribution of stars in the spheroidal component for our fiducial model if spheroid grow th through disk instability is 
suppr essed. Solid orange and green histograms sh ow observational measur ements bvlWvse fc Gilmord l|l995i . top left panel). iZoccali et al.l 
1I2OO3I . top right panel) , iNordstrom et al.l 1 120041 . bottom left panel) , and IZoccali et al.l 1120081 . bottom right p anel). The dashed orange 
histogram in the top left panel has been obtained converting the [Fe/H] scal e of the original distribution by IWvse fc Gilmord into an 
[O/H] scale by using the observed [0/H]-[Fe/H] relation for thin disk stars bv lBensbv. Feltzing fc LundstromI ||200J). 



6 THE STELLAR HALO 

As mentioned in Sec.[l] the stellar halo of our Galaxy has a 
mass that is less than about 1 per cent of the total stel- 
lar mass (i.e. about IO^'Mq). Most of its stars are very 
old (probably older than 12 Gyr) and metal-poor, with en- 
hanced values of the elemental abundance ratio [a/Fe] . The 
currently accepted scenario is that the stellar halo formed 
(at least in part) from stars stripped from satellite galax- 
ies that were accreted by the Galaxy over its lifetime (the 



ISearle fc ZinnI scenario discussed in Sec. [T]). In this section, 
we study the formation and structural properties of this 
galactic component, in the framework of the model discussed 
above. 



Our working hypothesis is that the stellar halo built up 
from the cores of the satellite galaxies that merged with the 
Milky Way (or better its main progenitor - see below) over 
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Figure 8. Full merger tree of our model Milky Way galaxy from the lowest resolution simulation used in this study (GAO). The area of 
the symbols scale with the number of particles of the associated halo, and their colour indicates different galaxy types: blue for galaxies 
sitting at the centre of the main halo (those on which cooling is allowed), green for satellites galaxies associated to distinct dark matter 
substructures, and red for 'orphan' galaxies (see text). The black circles mark the galaxies accreted onto the main branch of the model 
Galaxy, traced back until the last time they are central galaxies. 



its lifetimelj In order to identify the stars that end up in the 
stellar halo, we construct the full merger tree of our model 
Milky Way galaxy, and identify those galaxies that merge 
onto the 'main branch', i.e. the branch that is obtained by 
connecting the system to its most massive progenitor at eac h 
node of the tree (see Sec. 4.1 of iDe Lucia fc Blaizotll2007l ). 
We then trace back each of these galaxies until they are 
for the last time central galaxies, and save the information 
(identification numbers and positions) of a fixed fraction of 
the most bound particles of their parent haloes at this time. 
We refer to these as 'star particles', and 'tag' them with the 
stellar metallicity of the galaxies residing at their centres. 

The procedure outlined above is illustrated in Fig. [8] 
which shows the full merger tree of our model Milky- Way 
galaxy from the lowest resolution simulation used in our 
study (GAO). The Milky- Way galaxy is shown at the top 
of the plot, and all its progenitors and their histories are 
shown downward going back in time. In this figure, the area 
of the symbols scales with the number of particles in the 
associated halo, while different colours are used for differ- 
ent 'types' of galaxies. Central galaxies (the only ones on 
which cooling is allowed) are shown in blue. Green symbols 
indicate galaxies associated to a distinct dark matter sub- 
structure, and red symbols indicate 'orphan' galaxies, i.e. 
galaxies whose parent halo has been reduced below the res- 
olution limit of the simulation. For visualisation purposes, 

^ Note that some fraction of these galaxies may also contribute to 
the bulge component, particularly those with the highest density. 



we have kept constant the mass of the dark matter sub- 
structure when a galaxy becomes orphan. The black circles 
in Fig.[8]mark the galaxies accreted onto the main branch of 
the model Galaxy (the leftmost branch in the figure), traced 
back until the last time they are central galaxies. In the fol- 
lowing, we refer to this time as the time of accretion. Fig. [8] 
shows that the star particles that make up the stellar halo 
in GAO are old: the majority of the satellites building up 
the stellar halo were 'accreted' at lookback times larger than 
~ 9 Gyr (although about half of them 'merge' onto the main 
branch much later). Only one (relatively massive) satellite 
was accreted at lookback time ~ 8 Gyr, and it merges onto 
the main branch only at lookback time ~ 1 Gyr. For the 
simulations GA2 and GAS, all accretions occur at lookback 
times larger than ~ 11 Gyr. As we tag our star particles 
using the stellar populations of the central galaxies at the 
time of accretion, all stars in the stellar halo of these model 
galaxies have ages larger than 11 Gyr. 

Fig. [8] illustrates that once dark matter haloes are 
accreted onto a larger system, they survive as distinct 
dark matter substru ctures for a relatively short time (e.g 
IPe Lucia et al.ll2004l ). but the galaxies residing at their cen- 
tre merge onto the main progenitor of the Milky Way galaxy 
much later. As discussed in Sec. [31 our model galaxies are 
not affected by the tidal stripping and truncation that ef- 
ficiently reduces the mass of the dark matter substructure s 
IJGhJgna et al.ll2000l : IPe Lucia et"alll2004l : iGao et aLlliooJ ). 
If this effect is important, the survival time-scales of orphan 
galaxies, as well as the stellar mass of survived satellites are 
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Figure 9. Stellar masses (top panel) and metallicities (middle 
panel) for all galaxies accreted onto the main branch for the sim- 
ulation GAS, as a function of the lookback time of the galaxy's 
merger. The solid red lines in the top and middle panels show the 
evolution of the stellar mass and metallicity in the main progen- 
itor of the Milky Way. The bottom panel shows the number of 
particles associated to the dark matter halo at the time of accre- 
tion, again as a function of the merging time of the galaxy that 
is located at its centre. Red symbols indicate objects associated 
to substructures with more than 500 particles. Red open symbols 
correspond to red filled circles but are plotted as a function of the 
time of accretion. 



likely to be overestimated. Work is underway to calibrate 
stellar tidal stripping using hydrodynamic simulations. 

Fig. 1^ shows the stellar masses (top panel) and metal- 
licities (middle panel) of all galaxies accreted onto the main 
branch for the highest resolution simulation used in our 
study (GAS), as a function of the lookback time of the 
galaxy's merger. Note that stellar masses and metallicities 
correspond to those at the time of accretion. The solid red 
lines in these panels show the evolution of the stellar mass 
and of the stellar metallicity in the main progenitor of the 
Milky Way galaxy (as in Fig.jS)). The bottom panel of Fig. [9] 
shows the number of particles associated to the dark matter 
haloes before accretion, again as a function of the merging 
time of the galaxies that reside at their centre. Red sym- 



bols indicate objects belonging to haloes with more than 
500 bound particles before accretion. Note that for these 
objects we have indicated both the time of accretion (open 
symbols) and the merger time (solid symbols; given by the 
dynamical friction timescale as described above). The figure 
shows that most of the galaxies that merge onto the main 
branch have stellar masses and metallicities that are much 
smaller than the current mass of the main progenitor over 
most of the galcixy's life-time. Most of the accreted galax- 
ies lie in quite small haloes and only a handful of them are 
attached to relatively more massive systems. These are the 
galaxies that contribute most to the build-up of the stellar 
halo. The red symbols in the bottom panel of Fig. |9] show 
that most haloes with more than ^ 500 particles were all 
disrupted more than ~ 6 Gyr ago. These haloes contain 
the few galaxies with stellar mass larger than 10® Mq which 
merge onto the main branch over the galaxy's life-time (top 
panel). The stellar metallicities are generally relatively low, 
with a median value of ~ 0.3 Zq, with larger values asso- 
ciated to larger galaxies (see below). The results illustrated 
in Fig. 1^ are in good agreement with those by iFont et al.l 
(|2006r ) who find that one or a few more satellites in the range 
10* — 10^° M0 can make up 50-80 per cent of the stellar halo 
and most of them are accreted early on (taccr > 9 Gyr) . 



Fig. llOl shows the projected distribution of the star par- 
ticles that end up in the stellar halo of the highest resolution 
simulation used in our study (GAS). For this figure, and for 
the results presented in the following, we use 10 per cent of 
the most bound particles of every accreted halo, but sim- 
ilar results are obtained if 5 per cent of the particles are 
selected. We obtain 4188 star particles for the simulation 
GA2, and S1716 star particles for the simulation GAS. The 
star particles in Fig. [10] are colour coded as a function of 
their metallicity, as indicated in the top-left panel, and their 
spatial distribution is shown at the same redshifts used for 
Fig. [l] The star particles that end up in the stellar halo 
extend over a projected region of ~ 1 Mpc^ comoving at 
redshift ~ 10. At redshift ~ 1, the star particles are already 
assembled in a single relatively elongated component which 
becomes progressively more spherical with decreasing red- 
shift (see below). We note that the the 'stellar halo' of the 
simulation GAS is dominated by particles contributed by a 
single object that merged ~ 7 Gyr ago (and whose dark halo 
was accreted 11 Gyr ago) and with metallicity 0.5S Zq (see 
Fig. [9]). This single galaxy contributes 19921 star particles 
(i.e. about 60 per cent of all the star particles in the stellar 
halo) . 



The lower panels of Fig. [10] seems to suggest a slight 
concentration of the most metal-rich stars, but no clear cor- 
relation between metallicity and distance, with low metal- 
licity and high metallicity stars distributed all at various 
distances. It is interesting that in some cases there is a clear 
difference between the mean metallicity of stars in a given 
satellite (pockets of stars of different colours) - which is 
the basis of the 'chemical tagging ' argument discussed in 
[Freeman fc Bland-HawthornI (|2002l ). The resolution of the 
simulations used in this study is, however, too low to see any 
spatially coherent stellar streams in the halo at the present 
day. 
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Figure 10. Projected distribution of the star particles tliat end up in the stellar halo of the highest resolution simulation used in our 

study (GAS), at the same redshifts as in Fig. [T] and colour-coded as a function of their metallicity as indicated in the top left panel. If 

we assume a [a/Fe ]~ +0.4 dex, the highest metallicity bin could be translated into —0.7 < [Fe/H], while the lowest to [Fe/H]< —1.6 

jSalaris et al.l ll993^. As for Fig. [T] the box at 2 = is centred on the most bound particle of the GAS halo, while those at higher redshifts 

are centred on the most bound pa rticle of the main progenitor at the corresponding redshift. 
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Figure 11. Mean (filled circles) and median (empty circles) 
metallicity of star particles as a function of the distance from 
the most bound particle in the Milky Way halo for the simula- 
tions GA2 (blue) and GAS (red). Dashed lines correspond to the 
15th and 85th percentiles of the distributions. 



The absence of a clear metallicity gradientlj is illustrated 
more explicitly in Fig. [11] which shows the metallicity of the 
star particles as a function of the distance from the most 
bound particle in the Milky Way halo for the simulations 
GA2 (blue) and GAS (red). Filled and empty circles con- 
nected by solid lines indicate the mean and the median of 
the distributions respectively. Dashed lines correspond to 
the 15th and 85th percentiles. In both simulations, the mean 
metallicity decreases from Log[Z/ZQ] ~ —0.4 at the centre 
to ~ —0.8 at a distance of ~ 40 kpc. The median and upper 
85th percentile of both distributions are approximately flat 
around ~ —0.5. The lower 15th percentile exhibits a cer- 
tain decline with increasing distance from the centre, sug- 
gesting that the inner region is largely dominated by high- 
metallicity stars while the contribution from lower metallic- 
ity stars becomes more important moving to the outer re- 
gions. We recall, however, that both distributions are dom- 
inated in number by star particles associated to one or a 
few accreted galaxies with relatively high metallicity (hence 
the flat behaviour of the median and upper percentile): 58 
per cent of the steUar halo particles in GA2 have metallicity 
^ —0.4, and this fraction rises to 67 per cent for the simu- 
lation GAS. As discussed previously, the metallicity of our 
stellar halo is higher than what is known for the Galactic 
halo near the Sun. This could be reflecting that the metal- 
licity of the accreted galaxies may be too high (see below) . 

From the observational viewpoint, it is interesting to 
ask where the most metal poor stars of the stellar halo are 
located. We address this question in Fig. [12] which shows 
the projected density profiles (solid black lines) of the star 
particles for the simulations GA2 (top panel) and GAS (bot- 



We have verified that an age gradient is also absent. 
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Figure 12. Projected density profile of the stellar halo (solid 
black lines) and of the dark matter halo (dashed black lines) for 
the simulations GA2 (top panel) and GAS (bottom panel). The 
solid green and orange lines show the projected density profiles 
for star particles with metallicity smaller and larger than 0.4 Zq 
respectively. The blue and red dashed lines show Sersic fits to the 
profiles of star particles in the two metallicity bins. 



tom panel) obtained by stacking the three projections on 
the xy, xz, and yz plane. The dashed black lines in Fig. [12] 
show the projected density profile of th e dark matter halqj. 
In ag reement with previous findings (jBullock fc Johnston! 
[2OO5I), we find that the profile of the stellar halo is steeper 
and more centrally concentrated than the dark matter pro- 
file. The half-mass radius of the stellar halo in GAS is located 



^ We use only the particles in the 'main halo', i.e. the self-bound 
part of the friends-of-friends halo. 
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Table 3. Results of the Sersic fit (see text) to the projected 
profiles of star particles with metallicity larger and smaller than 
0.4 Zq (orange and green lines in Fig. I12I I. 



Figure 13. Stellar metallicity as a function of stellar mass for the 
galaxies contributing to the stellar halo in the simulation GAS. 



at ~ 4.3 kpc, i.e. well within the Solar c ircle, as estimated 
for the Milky Way (JFrenk fc Whitdll982l ). The solid orange 
and green lines show the projected profiles of star particles 
with metallicity larger and smaller than 0.4 Zq respectively. 
Fie:. ll2l shows that star particles with metallicity larger than 
0.4 Zq are more centrally concentrated than star particles of 
lower abundances. This implies that the probability of ob- 
serving low-metallicity stars increases if one looks at larger 
distances from the Galactic centre (> 10 — 20 kpc) where the 
contribution from the inner more metal-rich star particles is 
less dominant. 

Interestingly, our results appear to be in qualitativ e 
agreement with the measurements bv lCarollo et al.l ((20071). 
These authors have analysed a large sample of over 20,000 
stellar spectra from SDSS, and demonstrated that the stellar 
halo of the Milky Way can be considered as the superposi- 
tion of two components which are spatially, kinematically, 
and chemically distinct. One component, which they call 
'inner halo', dominates the population of stars found at dis- 
tances up to 10-15 kpc from the Galactic centre and peaks 
at higher metallicity than the 'outer halo', which dominates 
in regions beyond 15-20 kpc. 

In order to understand this 'duality' we turn to our sim- 
ulations. From the middle panel of Fig.|5]we see that stellar 
metallicity does not correlate with the time of the merger of 
the galaxies that contribute to the stellar halo. Rather, the 
'duality' of our model stellar halo originates from a correla- 
tion between the stellar metallicity and the mass in stars of 
the accreted galaxies, as shown in Fig. 1131 Since the most 
massive galaxies can decay through dynamical friction to the 
inner regions of the host halo, this is where higher metal- 
licity stars will be found preferentially. Note that the mass- 
metallicity relation in our simulations appears to be offset 
from that observed in the sense that the accreted galaxies 
are too metal-rich, explaining why our stellar halo also has 
too high a metallicity. This might be due to a too efficient 
mixing of metals with the intergalactic medium (we recall 
that we are using a 100 per cent mixing efficiency), and/or 
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to inefficient feedback in the building blocks of the stellar 
halo. We will analyse this in further detail in a future paper 
where we will study the satellite population of our model 
Milky Way. 

A large fraction of galaxies accreted onto the main 
branch contribute to the star particles in the inner halqf], 
i.e r < 10 kpc, but only a few systems contribute signifi- 
cantly to the population of star particles in this region. For 
both simulations, about 65 per cent of the star particles in 
the inner region comes from one single system that merged 
onto the main branch ~ 8 Gyr ago for the simulation GA2, 
and ~ 7 Gyr ago for the simulation GA3. Another ~ 20 per 
cent of the star particles in the inner 10 kpc are contributed 
in approximately equal fractions from other three galaxies 
that merged between ~ 8 and ~ 10.5 Gyr ago for the simu- 
lation GA2, and a bit earlier (between ~ 10 and ~ 11.5 Gyr 
ago) for the simulation GAS. 

The blue and red dashed lines in Fig. [T2] show Sersic 
fits to the star particles with metallicity larger and smaller 
than 0.4 Zq respectively (corresponding to the orange and 
green lines). The fitted function is given by: 

E(i?) = Eetf exp [-6„[(i?/i?eff )<'/"' - 1]] 

where 7?eff is the radius containing half the light and EoH 
is the surface brightness at that radius. The results of the 
fit are given in Table O For both simulations, the profile 
of star particles with metallicity lower than 0.4 Zq is well 
fit by a Sersic profile with index n ~ 3 and half-light ra- 
dius ~ 4 kpc. For the star particles with larger metallicity 
the characteristic radius is comparable, and the Sersic index 
is smaller (the distribution falls off more steeply) . The log- 
arithmic slope of the corresponding density profile at R^g 
is 7 ~ —3.3 for the metal-rich component and 7 ~ —3.1 
for the metal-poor stars, in very good agreement with those 
measured for the stellar halo of the Milky Way at a similar 
galactocentric distance. 

In order to characterise the three-dimensional shape of 
the stellar halo, we have assumed that its isodensity surfaces 
can be approximated by triaxial ellipsoids of the form: 



X? 



+ 



fe2 



+ 



xi 



— 1 with a ^ b ^ c 



where a, b, and c are the lengths of the three axes and X^ 
(with a — 1,2,3) is the coordinate with respect to the a 

** 70 and 95 per cent of the accreted galaxies contribute at least 
one particle to this region of the stellar halo for the simulations 
GA2 and GAS respectively. 
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Table 4. Axial ratios for the stellar halo and for the dark matter 
halo of the Milky Way galaxy for the simulations GA2 and GAS 
used in this study. 



stellar halo 


c/a 


c/b 


GA2 
GAS 


0.51 
0.45 


0.74 
0.89 


dark matter halo 


c/a 


c/b 


GA2 
GAS 


0.70 
0.72 


0.80 
0.85 



axis and relative to the position of the most bound particle 
of the Milky Way halo. The directions and lengths of the 
major axes can then be computed by finding the eigenvectors 
and eigenvalues of the matrix Ma/s: 



M, 



al3 



7 , X^X/j 



where a,f3 — 1, 2, 3. In Tabled we list the resulting axial ra- 
tios c/a and c/b for the stellar halo and for the dark matter 
halo, for the simulations GA2 and GAS. The principal axes 
of the stellar halo and of the dark matter halo are plotted 
as thin and thick solid lines in Fig. 1141 The top panel cor- 
responds to the simulation GA2, while the bottom panel is 
for the simulations GAS. In both panels, orange and green 
points are star particles with metallicity larger and smaller 
than 0.4 Zq 

The value of the short-to- long (c/a) axis ratio measured 
for the dark matter halo is in good agreemen t with p r evious 
numerical studies in a ACDM cosmology bv iBuUockl (120021 ) 
who finds typical axis ratios to be c/a ~ 0.6 — 0.8 for Milky 
Way size haloes. The stellar halo has a short-to-long axis ra- 
tio ~ .45 — 0.5, slightly lower than observat ional measure- 
ments (|Chiba fc BeerslbOOd : iBeU et ahllJOOTl , and references 
therein). Interestingly, for both simulations the minor axis of 
the stellar halo is quite well aligned with the the minor axis 
of the dark matter halo: the angle between these two axes 
is only ~ 10 deg for the simulation GA2 and ~ SO deg for 
the simulation GAS. For the highest resolution simulation 
(GAS), also the major and intermediate axes of the stellar 
halo are well aligned with the corresponding axes of the dark 
matter halo (the angular distance is ~ 9 and ~ SO deg for 
the other two axes). For the simulation GA2, the angular 
distance between the major and intermediate axes of the 
stellar halo and the corresponding axes of the main halo is 
~ 75 deg. 



7 DISCUSSION AND CONCLUSIONS 

During the last decade, a number of observational tests of 
the standard cosmological model have ushered in a new 
era of 'precision cosmology'. While the basic theoretical 
paradigm for structure formation is now well established, our 
understanding of the physical processes governing galaxy 
formation and evolution is far from complete. Modern ad- 
vances in ground- and space-based observational capabilities 
have enabled us to study the resolved stellar populations of 
Local Group galaxies (and beyond) in unprecedented detail. 
For our own Galaxy and for a number of the brightest mem- 





Figure 14. Three-dimensional distribution of the star particles 
from the simulations GA2 (top panel) and GAS (bottom panel). 
Green symbols refer to particles with Z/Zq < 0.4, while orange 
symbols refer to star particles with Z/Zq J^ 0.4. Red, blue and 
green lines represent the principal axes of the stellar halo (thin 
lines) and of the dark matter halo (thick lines). Red is used for 
the long axis, blue for the intermediate, and green for the short 



bers of the Local Group, a wealth of observational data is 
now available about the ages and chemical abundances of 
their stars. Much more data will become available over the 
next decade, providing important tests for current models 
of galaxy formation and evolution. 

In this paper, we have discussed the formation of the 
Milky Way and of its stellar halo in the context of a hybrid 
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cosmological approach that combines high-resolution simu- 
lations of a 'Milky- Way' halo with semi-analytic methods. 

Our approach is s imilar to that adopted by 
iBullock fc JohnstonI (120051). but the two m ethods differ in 
a number of details. iBullock fc JohnstonI have used mass 
accretion histories of Milky Way-size galaxies using the ex- 
tended Press-Schechter formalism. For each accretion event 
in the analytic merger trees, they have run A''-body simula- 
tions following the dynamical evolution of the accreted dark 
matter satellite in an analytic parent galaxy-|-host halo po- 
tential. The stellar distribution and properties of each satel- 
lite are constructed using a variable mgiss -to-light ratio to 
each dark matter particle IjFont et al.ll2006l ). and the chemi- 
cal evolution of the satellites is modelled taking into account 
the e nrichment from both Type II and Type la supernovae 

I Rob ert son ct al. 2005). In our work, we have used a set of 
increasing resolution A'^-body simulations of a Milky Way- 
like halo, combined with a semi-analytic method that al- 
lows us to model the stellar distribution of the galaxies self- 
consistently during the N-body simulation. While our ap- 
proach is based on a fully numerical simulation, the numer- 
ical resoluti on for each accreted satellite is lower than that 
achieved bv lBuUock fc JohnstonI . In addition, our chemical 
evolution model assumes an instantaneous recycling approx- 
imation which is appropriate for elements produced by Type 

II supernovae, but wrong for the iron-peak elements which 
are mainly produced by Type la supernovae. 

The galaxy formation model employed in this work has 
been studied in a number of previous papers, and it has 
been shown to successfully reproduce a number of observa- 
tional results for the global galaxy population in the local 
Universe and at higher redshifts. The physical properties of 
our model Milky- Way galaxies are in quite nice agreement 
with the observational results, and the predicted evolutions 
are very similar, over the entire numerical range covered 
by the simulations used in our work (from ~ 2 x 10* to 
~ 3 X 10^ Mq dark matter particle mass) . This is the first 
time to our knowledge that model convergence has been 
shown on a galaxy- by- galaxy basiqj. 

Our model Milky Way galaxy is a relatively young sys- 
tem (50 per cent of the the stars are assembled in a single 
object only at redshift ~ 0.8) with a very old spheroidal 
component (50 per cent of the stars in the spheroid are in 
place at redshift ~ 3), which is formed through a series of 
minor mergers and a few episodes of disk instability occur- 
ring early on during the galaxy's lifetime. All stars in the 
spheroidal component are old (> 11 Gyr) while the stars in 
the disk have a much larger spread in age, reflecting a pro- 
longed star formation activity which is in qualitative agree- 
ment with observational determinations. A detailed compar- 
ison between model and observed metallicity distributions 
is complicated by the use of the instantaneous recycling ap- 
proximation, which does not allow direct comparison with 
the observed iron distribution. Using the observed relation 
between oxygen and iron abundances, we obtain a metallic- 
ity distribution which is similar to the observed one for the 
disk component, but a spheroidal component which is less 



enriched than the observed Galactic bulge (although small 
changes in the parameters of the semi-analytic models are 
able to remove this discrepancy without signiflcantly affect- 
ing the remaining properties of the simulated Galaxy, see 
e.gFig.HJ. 

Assuming that the stellar halo builds up from the cores 
of the satellite galaxies that merged with the Milky Way 
over its life-time, we are able to study the physical and 
structural properties of this Galactic component. Our model 
stellar halo is made up of very old stars (the majority of the 
stars formed earlier than 10 Gyr ago) with low metallic- 
ity (^ 0.5 Zq), although relatively high in comparison to 
the Galactic halo near the Sun. A few relatively massive 
(10* - 10^" Mq) satellites accreted early (> 9 Gyr) con- 
tribute the largest fraction of the star particles that end 
up in th e stel l ar ha lo, in agreement with previous results 
bv lFont et al.l (|2006| V There is no evidence of a metallicity 
gradient for halo stars, but we find evidence for a stronger 
concentration of higher metallicity stars. This implies that 
the probability of observing low-metallicity halo stars in- 
creases with distance from the Galactic centre (> 20 kpc). 
The 'duality' we find for our model stellar halo does not 
originate from a correlation between metallicity and accre- 
tion time, rather from the fact that the main contributors 
to the stellar halo (satellites with larger mass) have pref- 
erentially higher metallicity than lower mass systems. We 
find that the three-dimensional distribution of halo stars is 
well described by a triaxial ellipsoid with short-to-long axis 
~ 0.5 — 0.6, in agreement with observational measurements, 
and whose axes are well aligned to those of the parent dark 
matter halo. 

A more accurate treatment of chemical enrichment will 
allow us to carry out a more detailed comparison with ob- 
served chemical compositions, and to establish similarities 
and differences between present-day satellites and the Galac- 
tic building blocks. The numerical resolution of the simu- 
lations used in this work, however, is too low for studies 
of spatially and kinematically coherent stellar streams in 
the present day stellar halo. Higher resolution simulations 
are therefore needed for these studies. These are all much 
needed steps to interpret the outcome of large surveys such 
as SEGUE, RAVE and ultimately Gaia, with the goal of 
unrevealing the evolutionary history of our Galaxy. 
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